class: title-slide, middle, right <br> <br> <br> .line_space_15[ ## .text_80[Modelos estructurales marginales para el control de<br>sesgos en estudios observacionales con factores de<br>riesgo y exposición tiempo-dependientes] ] <br> .line_space_11[ <br> .text_70[[Código en: <svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg>](https://github.com/AGSCL/DSPUCH)] .text_110[Seminario Métodos de Investigación en Salud Pública] ] .bg-text[ 10 de julio, 2022 .text_100[Andrés González Santa Cruz] .text_50[gonzalez.santacruz.andres@gmail.com] [<svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg>](https://github.com/AGSCL) [<svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:green;overflow:visible;position:relative;"><path d="M294.75 188.19h-45.92V342h47.47c67.62 0 83.12-51.34 83.12-76.91 0-41.64-26.54-76.9-84.67-76.9zM256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zm-80.79 360.76h-29.84v-207.5h29.84zm-14.92-231.14a19.57 19.57 0 1 1 19.57-19.57 19.64 19.64 0 0 1-19.57 19.57zM300 369h-81V161.26h80.6c76.73 0 110.44 54.83 110.44 103.85C410 318.39 368.38 369 300 369z"/></svg>](https://orcid.org/0000-0002-5166-9121) ] <img src="data:image/png;base64,#./_style/Logo_nDP_monotono_vertical_en.png" width="15%" /> ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Mi nombre es andrés gonzález y en esta instancia presentaré sobre los modelos estructurales marginales - Decidí abordar un tema metodológico no tratado en clases, que les será útil para mi tesis y que podría ser útil para otros. - --- layout: true class: animated, fadeIn --- ## Índice <br> 1. Problemática 2. Causalidad 3. Ajuste confusión 4. Marginal Structural Models 5. Aplicación 6. Desafíos ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA INDICE** *#_#_#_#_#_#_#_#_#_#_ Primeras dos, repaso: 1- **Origen de mi interés y problemática**: factores de riesgo y exposición tiempo-dependientes: Ejemplo de arbol de decisiones y que al final obtenía puros efectos condicionales. 2- **Causalidad**: La tarea de la epidemiología, definicion de causalidad y de confusión. Distribución marginal vs. conjunta, etc. Datos observacionales. Supuestos: SUTVA; Sesgos.Supuestos a la base, amenazas, sesgos, y particularmente sesgos que ocurren para este diseño de datos. 3- **Ajuste confusión**: Las formas. En una de esas, mostrar gif y animaciones (https://youtu.be/j8J2L_g76c4?t=68 o https://twitter.com/nickchk/status/1068215492458905600 / https://github.com/NickCH-K/causalgraphs/blob/master/Animation%20of%20IV.R / https://nickchk.com/causalgraphs.html). Explicar las diferencias entre estandarización y restricción --> ir poniendo el código paso por paso 4- **MSM**, los tipos.Explicar por qué son **estructurales** (1) y **marginales**(2), explicar las situaciones en que se ocupan y son útiles. Incluso con datos aleatorizados. Sesgos que resuelven 5- **Aplicación**: Explicar estructura hipotética de los datos, aplicación de un modelo simple. 6- **Desafíos y aproximaciones más complejas** --- ## Problemática .panelset.sideways[ .panel[.panel-name[Estructura hipotética más cercana a lo real] .details-code[ ```r # Libraries library(ggplot2) set.seed(2125) # Create data data <- data.frame( y=abs(rpois(1:250,15)), y2=abs(rpois(1:250,15)) ) %>% dplyr::filter(y2>y, y>=8, y2<=22) %>% dplyr::mutate(Paciente=row_number()) %>% #filtrar tratamientos más largos que 3 años dplyr::mutate(diff_treat=y2-y) %>% dplyr::filter(diff_treat<=3) for (i in 1:nrow(data)){ data$y3[i]<-base::sample(x=seq(from=data$y2[i]+1,to=22),1) data$y3[i]<-ifelse(data$y3[i]<=data$y2[i],22,data$y3[i]) data$y3[i]<-ifelse(data$y3[i]>=23,22,data$y3[i]) data$y4[i]<-ifelse(!is.na(data$y3[i]),base::sample(x=seq(from=data$y3[i]+1,to=22),1),22) data$y4[i]<-ifelse(data$y4[i]<=data$y3[i],22,data$y4[i]) data$y4[i]<-ifelse(data$y4[i]>=23,22,data$y4[i]) } set.seed(2125) pac_aleatorio1<-sample(1:max(data$Paciente),40) pac_aleatorio2<-setdiff(sample(1:max(data$Paciente),40), pac_aleatorio1) set.seed(2125) pac_aleatorio3<-unique(data$Paciente)[sample(1:length(data$Paciente),15)] # Horizontal version, antes era 1985 end_plot<-20 fig_trans<-ggplot(data) ###>=10 #datos de tratamiento en periodo de seguimiento, posterior al 2010 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2<=end_plot, Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=y, yend=y2), color="#EF9D2F", alpha=.6, size=.8)) fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2<=end_plot, !Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=y, yend=y2), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entera hasta el 2020 fig_trans<-try(fig_trans+ geom_segment(data=dplyr::filter(data,y>=10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=y, yend=end_plot), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=end_plot, yend=y2), color="#21177A", alpha=.6, linetype="dotted", size=.8)) ###<10 #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10), aes(x=Paciente, xend=Paciente, y=y, yend=10), color="#21177A", alpha=.6, linetype="dotted", size=.8)) #datos de tratamiento en periodo de seguimiento, menos de 2010 pero sólo y1 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>=10,y2<=end_plot, Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=10, yend=y2), color="#EF9D2F", alpha=.6, size=.8)) #fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>=10,y2<=end_plot, !Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=10, yend=y2), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entera hasta el 2020 # fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=10, yend=end_plot), color="#21177A", alpha=.6, size=.8)) # CONTACTO: VICTIMARIO: debo saber si y3 está dentro o no del seguimiento ###>=10 #datos de justicia en periodo de seguimiento, linea solida fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4<=end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=y4), color="#21177A", alpha=.6, size=.8)) #datos de justicia en periodo de seguimiento, exceden el seguiemiento pero parten en él. linea solida hasta el 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=end_plot), color="#21177A", alpha=.6, size=.8)) #si estoy dentro de los datos en y, peroen end plot no, tengo q hacer el interlineado afuera fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=end_plot, yend=y4), color="#21177A", alpha=.6,linetype="dotted", size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3>end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=y4), color="#21177A", alpha=.6, linetype="dotted", size=.8)) fig_trans_final<- #Puntos fig_trans+ #Tratamientos completados geom_point(data=dplyr::filter(data,Paciente%in% sample(c(pac_aleatorio1,pac_aleatorio2),15), Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y2), color="#EF9D2F", size=3, alpha=.6) + #geom_point(data=dplyr::filter(data,Paciente%in% pac_aleatorio1, Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y4),color="#EF9D2F",size=3,alpha=.6)+ geom_point(data=dplyr::filter(data,Paciente%in% sample(c(pac_aleatorio1,pac_aleatorio2),15), !Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y2), color="#21177A", size=3, alpha=.6) + geom_point(data=dplyr::filter(data,Paciente%in% pac_aleatorio1, !Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y4),color="#21177A",size=3,alpha=.6)+ theme_light() + coord_flip() + #fechas donde yo tomo gente annotate("rect", xmin=-Inf, xmax=Inf, ymin=10, ymax=10.5, alpha = .2, fill="blue")+ annotate("rect", xmin=-Inf, xmax=Inf, ymin=19.5, ymax=20, alpha = .2, fill="blue")+ theme( panel.grid.major.y = element_blank(), panel.border = element_blank(), axis.ticks.y = element_blank() )+ scale_x_continuous(breaks=seq(1,max(data$Paciente),by=10))+ scale_y_continuous(breaks=seq(min(data$y),max(data$y4),by=2), labels=seq(min(data$y),max(data$y4),by=2)+2000)+ labs(y="Tiempo de seguimiento (en trimestres)", x="Individuos (ID)")+ #labs(y="Follow-up (in quarters)",x="Individual(ID)", caption="Note. Dot= Complete treatment;\nSquare= Contact w/ justice system (imputed);\nTriangle= Contact w/ justice system (victim);\nBlue line= Time in treatment;\nOrange line= Time in conctact w/justice system;\nShared area=Follow-up window")+ theme(plot.caption = element_text(hjust = 0, face= "italic"))+ labs(caption="Punto= Completa treatmento;\nAzul= Tratamiento ambulatorio;\nAmarillo= Tratamiento residencial; \nLínea= Tiempo en tratamiento;\nArea sombreada=Inicio y término seguimiento")+ theme( panel.background = element_rect(fill = "transparent", colour = NA_character_), # necessary to avoid drawing panel outline #panel.grid.major = element_blank(), # get rid of major grid #panel.grid.minor = element_blank(), # get rid of minor grid plot.background = element_rect(fill = "transparent", colour = NA_character_), # necessary to avoid drawing plot outline legend.background = element_rect(fill = "transparent"), legend.box.background = element_rect(fill = "transparent"), legend.key = element_rect(fill = "transparent")) invisible("Paciente 49: 2011 al 2014 cursa un primer tratamiento que abandona, luego del 2016 al 2018 completó un tratamiento") fig_trans_final ``` <img src="data:image/png;base64,#seminario_files/figure-html/p1-est-1.png" width="70%" height="70%" style="display: block; margin: auto;" /> ] ] .panel[.panel-name[Primera aproximación] Objetivo: **Estimar el efecto de la modalidad de tratamiento por trastornos por uso de sustancias (ambulatorios vs. residenciales) a la base en la probabilidad de experimentar subsiguientes readmisiones a TUS** Análisis propuestos en esa oportunidad: .superbigimage[ .text_50[ `$$readmisión\sim A_{1 (Modalidad\,base)}+X_{2 (Completa\,tto.\,base)}+\epsilon$$` `$$2^{da}\,readmisión\sim A_{1(Modalidad\,base)}+L_{2 (Completa\,tto.\,base)}+L_{3 (Completa\,2^{do}\,tto.)}+L_{4 (Días\,desde\,1^{er}\,ingreso)}+ \epsilon$$` `$$3^{era}\,readmisión\sim A_{1(Modalidad\,base)}+ L_{2 (Completa\,tto.\,base)}+ L_{3 (Completa\,2^{do}\,tto.)}+ L_{4 (Completa\,3^{er}\,tto.)}+ L_{5 (Días\,desde\,1^{er}\,ingreso)}+\epsilon$$` `$$4^{ta}\,readmisión\sim A_{1 (Modalidad\,base)}+ L_{2 (Completa\,tto.\,base)}+ L_{3 (Completa\,2^{do}\,tto.)}+ L_{4 (Completa\,3^{er}\,tto.)}+ L_{5 (Completa\,4^{to}\,tto.)}+ L_{6 (Días\,desde\,1^{er}\,ingreso)}+\epsilon$$` ] ] .center[**¿Cómo se reflejaría en un diagrama?**] ] .panel[.panel-name[Diagrama Causal Simplificado] .details-code[ ```r dag23 <- dagitty('dag { bb="0,0,1,1" "A0 (Modalidad tto. basal)" [exposure,pos="0.184,0.608"] "A1 (Modalidad tto. 2)" [pos="0.451,0.605"] "A2 (Modalidad tto. 3)" [pos="0.734,0.605"] "L0 (Vector car. individuales)- Pareamiento" [adjusted,pos="0.121,0.397"] "L1 (Días previos a la readmisión 2)" [adjusted,pos="0.507,0.341"] "L2 (Días previos a la readmisión 3)" [adjusted,pos="0.691,0.365"] "LM0(Completa tto. basal)" [adjusted,pos="0.269,0.707"] "LM1 (completa tto. 2)" [adjusted,pos="0.566,0.714"] "LM2 (completa tto. 3)" [adjusted,pos="0.838,0.722"] "U0 (ej., factores genéticos)" [latent,pos="0.124,0.273"] "U1 Car. indivudales (ej., patrón de consumo)" [latent,pos="0.444,0.829"] "U2 Car. individuales T2" [latent,pos="0.733,0.814"] "Y1(Readmisión)" [outcome,pos="0.355,0.604"] "Y2(2da Readmisión)" [outcome,pos="0.626,0.605"] "Y3(4ta Readmisión)" [outcome,pos="0.960,0.607"] "A0 (Modalidad tto. basal)" -> "LM0(Completa tto. basal)" "A0 (Modalidad tto. basal)" -> "Y3(4ta Readmisión)" "A1 (Modalidad tto. 2)" -> "LM1 (completa tto. 2)" "A2 (Modalidad tto. 3)" -> "LM2 (completa tto. 3)" "L0 (Vector car. individuales)- Pareamiento" -> "A0 (Modalidad tto. basal)" "L0 (Vector car. individuales)- Pareamiento" -> "Y1(Readmisión)" "L1 (Días previos a la readmisión 2)" -> "Y2(2da Readmisión)" "L2 (Días previos a la readmisión 3)" -> "Y3(4ta Readmisión)" "LM0(Completa tto. basal)" -> "U1 Car. indivudales (ej., patrón de consumo)" "LM0(Completa tto. basal)" -> "Y1(Readmisión)" "LM0(Completa tto. basal)" -> "Y2(2da Readmisión)" "LM0(Completa tto. basal)" -> "Y3(4ta Readmisión)" "LM1 (completa tto. 2)" -> "U2 Car. individuales T2" "LM1 (completa tto. 2)" -> "Y2(2da Readmisión)" "LM1 (completa tto. 2)" -> "Y3(4ta Readmisión)" "LM2 (completa tto. 3)" -> "Y3(4ta Readmisión)" "U0 (ej., factores genéticos)" -> "L0 (Vector car. individuales)- Pareamiento" "U0 (ej., factores genéticos)" -> "Y1(Readmisión)" "U0 (ej., factores genéticos)" -> "Y2(2da Readmisión)" "U0 (ej., factores genéticos)" -> "Y3(4ta Readmisión)" "U1 Car. indivudales (ej., patrón de consumo)" -> "A1 (Modalidad tto. 2)" "U1 Car. indivudales (ej., patrón de consumo)" -> "LM1 (completa tto. 2)" "U2 Car. individuales T2" -> "A2 (Modalidad tto. 3)" "U2 Car. individuales T2" -> "LM2 (completa tto. 3)" "Y1(Readmisión)" -> "A1 (Modalidad tto. 2)" "Y1(Readmisión)" -> "U1 Car. indivudales (ej., patrón de consumo)" "Y2(2da Readmisión)" -> "A2 (Modalidad tto. 3)" "Y2(2da Readmisión)" -> "U2 Car. individuales T2" }') tidy_dag23 <- tidy_dagitty(dag23) %>% dplyr::mutate(label=dplyr::case_when(grepl("A0",as.character(name))~"A0", grepl("A1",as.character(name))~"A1", grepl("A2",as.character(name))~"A2", grepl("L0",as.character(name))~"L0", grepl("L1",as.character(name))~"L1", grepl("L2",as.character(name))~"L2", grepl("LM0",as.character(name))~"LM0", grepl("LM1",as.character(name))~"LM1", grepl("LM2",as.character(name))~"LM2", grepl("U0",as.character(name))~"U0", grepl("U1",as.character(name))~"U1", grepl("U2",as.character(name))~"U2", grepl("Y1",as.character(name))~"Y1", grepl("Y2",as.character(name))~"Y2", grepl("Y3",as.character(name))~"Y3", T~as.character(name))) %>% dplyr::mutate(label2=dplyr::case_when(grepl("U",name)~"latent",grepl("L|LM",name)~"adj",grepl("LM",name)~"white", T~"black")) %>% dplyr::mutate(adjusted=factor(dplyr::case_when(grepl("L|LM",name)~"adjusted",T~"unadjusted"))) edge_function <- ggdag:::edge_type_switch("link_arc") dag23_plot<- ggdag:::if_not_tidy_daggity(tidy_dag23) %>% ggdag:::node_status() %>% ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend, yend = yend, color = status, shape=factor(adjusted)))+ #edge_function()+ scale_adjusted()+ ggdag:::breaks(c("exposure", "outcome","latent"))+ geom_dag_edges_arc(curvature = c(0,5,rep(0,26)))+ ggdag:::geom_dag_point(size = 16)+ ggdag:::geom_dag_label_repel(ggplot2::aes_string(label = "label", fill = "status"), size = 4.88, col = "white", show.legend = FALSE)+ theme_dag()+ scale_shape_manual(values = c(15, 16), name="Ajustado", labels=c("Sí", "No"))+ scale_fill_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado","No observado"), limits = c('exposure', 'outcome','latent'))+ scale_color_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado","No observado"), limits = c('exposure', 'outcome','latent'))+#E6E6E6 guides(linetype="none", edge_alpha="none", shape="none")+ guides(color=guide_legend(override.aes = list(arrow = NULL)))+#,guide_colourbar(order = 1) theme(plot.caption = element_text(hjust = 0))+ theme(legend.position = "bottom", aspect.ratio=6/10)+ labs(caption="Nota. Ak= Modalidad (Residencial/Ambulatoria); LMk= Completa tratamiento;\nL0= Confusores, características individuales y del centro a la base;Lk+= Días previos a la readmisión;\nYk= Resultadotratamiento; Uk=Car. no observadas para cada tiempo") dag23_plot ``` <img src="data:image/png;base64,#seminario_files/figure-html/p2-dag-1.png" height="450" style="display: block; margin: auto;" /> ] <div class="center"> ¿Qué ocurre con las modalidades en los tratamientos sucesivos?, ¿Confunden la asociación entre `\(A_0\)` ` e` `\(Y_3\)` `? `, ¿Sigue siendo causal en presencia de { `\(A_1\,, A_2\,,A_3\)` `}? </div> <!--- U puede tener relación con derivaciones no observadas posterior al tratamiento, que sería como un U1.5---> <!--- U puede tener relación con derivaciones no observadas posterior al tratamiento, que sería como un U1.5---> ] .panel[.panel-name[Árbol causal] <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-2-1.png" width="550" style="display: block; margin: auto;" /> ] .panel[.panel-name[Entonces...] - Varias puertas traseras sin cerrar, confusión residual, abandono/pérdida, sobreajuste, etc. - El diseño y la estrategia analítica hasta el momento no permite responder a una pregunta causal de manera adecuada (muchos contrastes condicionales a la trayectoria de tratamientos y del mediador) - Debiesen analizarse otras estrategias que capturen efectos longitudinales <br> .pull_r_30[ .text_60[ ***"In longitudinal studies with time-dependent confounding, identifying the structure allows us to detect situations in which stratification-based methods would adjust for confounding at the expense of introducing selection bias (p. 622)"*** `\(^{[1]}\)` ] ] <br> <br> <br> .down_center[ <div class="red"> ¿Se puede responder a la pregunta de si haber asistido a tratamiento residencial a la base ( `\(A_0\)` `) es beneficioso?, y si lo es, ¿qué régimen/estrategia es óptima o casi óptima? `\(^{[2]}\)` </div> ] ] ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ 1. Mi interés principalmente proviene de una investigación en la que estoy participando, en que analizo trayectorias de tratamientos de drogas del 2010-2019. En esta imagen sólo expongo una muestra con una estructura hipotética Cada fila es un paciente, una línea es un tratamiento, el espacio blanco entre ellos es el tiempo en que no están en tratamientos. Los puntos son completar el tratamiento y las líneas amarillas es cursar un tratamiento residencial. Por ejemplo, en este caso, el Paciente N°55 completa tto. residencial y experimenta el evento de interés (readmisión), en ambos completa el tratamiento. 2. El objetivo de ese estudio era estimar ese efecto mediante 4 regresiones secuenciales y estratificado por la presencia de cada evento (ej., se partió con 22.000 usuarios que sólo tenían una readmisión; y ya a la 4ta se estimó el efecto sólo con 200 usuarios). Como había una razonable sospecha de dependencia en tiempos anteriores, se generó un modelo semi-markov. Al margen de eso, hay varios problemas. 3. El modelo causal es difícil de representar para mí. Sobre todo porque a la base se generó un pareamiento. Hernan señala que no es posible representar en un diagrama un pareamiento (Al margen de eso, ya también es discutible su implicancia causal, ya que agrupa al 90% de los tratamientos residenciales, pero al 18% de los ambulatorios. Es un efecto en los tratados. Se puede observar que hay muchos colliders (en LM) y puertas traseras no cerradas ().The total effect cannot be estimated due to adjustment for an intermediate or a descendant of an intermediate.). 4. Robins 1986 en su momento estableció un gráfico bastante similar, en que L=0 y L=1 son los resultados del tratamiento (completa/no-completa), y los A son residencial/ambulatorio. Y= es readmisión  POR LO VISTO, SON EFECTOS ANIDADOS PRODUCTO DE LA COMBINACIÓN DE CADA CONFIGURACION (EJ., QUIEN PRESENTA UN REGIMEN EN QUE COMPLETA UN RESIDENCIAL A LA BASE, SEGUIDO DE COMPLETAR EL SEGUNDO TRATAMIENTO RESIDENCIAL,VA A TENER UN RIESGO DE READMISIÓN DISTINTO  DE AHÍ QUE EL EFECTO ES CONDICIONAL A ESTE REGIMEN DE EXPOSICION Y MEDIACION **Hasta el momento es encesario hacer modificaciones al diseño y la estrategia analítica** (por si ... ) - Que no haya Ajuste por variables intermedias o descendientes de una variable intermedia ("the mantra not to control for factors affected by exposure".- Overadjustment)Schisterman2009 - Control for `W` by seeing what `W` explains (perhaps with a regression) and taking it out - *crossover effects --- class: center, middle ## Objetivo Este proyecto apunta a servir como un material introductorio sobre modelos estructurales marginales para el control de sesgos en estudios observacionales con factores de riesgo y/o exposición tiempo-dependientes. -- - Se mostrará su estimación paso a paso y mediante un ejercicio reproducible. -- - Se utilizará una base de datos simulada para ilustrar el proceso. ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ Ejemplos de publicaciones de este tipo hay: - Tutoriales de herramientas analíticas innovadoras en "Psychological Methods" (APA) https://www.apa.org/pubs/journals/met/call-for-papers-tutorials - Bradburn, M., Clark, T., Love, S. et al. Survival Analysis Part I: Multivariate data analysis – an introduction to concepts and methods. British Journal of Cancer 89, 431–436 (2003). https://doi.org/10.1038/sj.bjc.6601119 - Zhang, Z., Abarda, A., Contractor, A., Wang, J., & Dayton, C. (2018). Exploring heterogeneity in clinical trials with latent class analysis. Annals Of Translational Medicine, 6(7), 9. doi:10.21037/atm.2018.01.24 --- ## Causalidad (1) .pull_left[ - Salud pública- epidemiología `\(^{[3]}\)` - Qué pasaría si... // Qué causa ... `\(^{[4; 5]}\)` <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-3-1.png" width="350" style="display: block; margin: auto auto auto 0;" /> - Asociación vs. Causalidad `\(^{[4; 6]}\)` - Inferencia `\(^{[6]}\)` - RCTs: gold-standard `\(^{[7]}\)` - Estudios observacionales .text_70[ ***"It’s much easier to get a result than it is to get an answer"*** `\(^{[8]}\)` ] ] .pull_right[ . <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-4-1.png" width="350" style="display: block; margin: auto 0 auto auto;" /> .line_space_05[ .text_60[ .footnote[ Riseberg,E., Melamed, R., James, K., Alderete, T. & Corlin, L. (2021). Development and application of an evidence-based directed acyclic graph to evaluate the associations between metal mixtures and cardiometabolic outcomes. doi: https://doi.org/10.1101/2021.03.05.21252993 ] ] ] ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** (1) *#_#_#_#_#_#_#_#_#_#_ 2- **Causalidad**: La tarea de la epidemiología, definición de causalidad y de confusión. Distribución marginal vs. conjunta, etc. Datos observacionales. Supuestos: SUTVA; Sesgos.Supuestos a la base, amenazas, sesgos, y particularmente sesgos que ocurren para este diseño de datos. - Muchas preguntas en salud pública se remiten a cuestiones epidemiológicas. Desde determinantes (ej., vulnerabilidad, PP PP), pasando por factores de riesgo, procesos biológicos/mecanismos etiológicos, etc. - Una asociación no siempre es causal, en esta presentación me remitiré a preguntas que apuntan a la identificación o a cuantificar relaciones causales. o Hay asociación cuando el tener información sobre una variable me permite predecir mejor en promedio el resultado de Y, incluso si esa variable no causa Y. El mecanismo no lo sabemos. o Uno de los desafíos es poder determinar si por sí sola una variable X causa Y - Las preguntas epidemiológicas no se limitan a describir sino que algunas apuntan a un análisis causal (Kaufman Santiago 2020 Lunes). Al margen de la tarea valorable de la vigilancia, la predicción en el futuro, el “Mundo hipotético”, ¿qué pasaría si hacemos/aplicamos X? (efectos de causas)/ qué causa Y (causas de efectos)” (Gelman2013) - Nos restringiremos a al estudio de las intervenciones. [Tener presente los alcances del estudio de potential outcomes framework] - En un mundo complejo como el de hoy, es necesario observar los fenómenos de manera multivariable, razonamiento estadístico, el efecto de intervenir sobre una variable sobre otra. Cummiskey2020 - Ensayos clínicos controlados, buscan estimar el impacto de un tratamiento, patrón oro de evidencia, aleatorizados (NEFRO) - Incluso ellos tienen limitaciones: no son aplicables, pueden carecer de potencia estadística, pueden ser no éticos, por diseño; por eso estudios observacionales (NEFRO) - El control de la confusión es uno de los pasos más desafiantes en la conducción de estudios observacionales (MANSOURNIA) y a menudo puede ocurrir también en ensayos clínicos. - El marco de inferencia contrafactual (potential outcomes framework) (Neyman, 1923; Rubin, 1974), iniciado por Neyman y continuado por Rubin y Robins, busca comparar teóricamente qué ocurriría si todos los pacientes hubieran sido tratados frente a si ningún paciente hubiese sido tratado. - Una asociación puede definirse en términos de una distribución conjunta de variables observadas, mientras que un concepto causal no puede definirse solamente de la distribución. (Pearl 2003). o Asociaciones: correlaciones, regresiones, dependencia, independencia condicional, verosimilitud, colapsabilidad, PS, RR, OR, marginalización, condicionalización, controlando por. o Causal: Aleatorización, influencia, efecto, confusión, manteniendo constante, disturbance, correlación espúrea, estabilidad, IV, intervención, explicación, atribución o “different risk in two disjoint subsets of the population determined 12 Causal Inference by the subjects’ actual treatment value ( = 1 or  = 0), whereas causation is defined by a different risk in the entire population under two different The difference between association treatment values ( = 1 or  = 0)” (Hernan y Robins) - La cita se vincula con que, para Pearl, el dato de por sí no es suficiente para señalar que una relación es causal(Pearl, 2003) Y eso nos da paso a la siguiente diapositiva. o These considerations imply that the slogan “correlation does not imply causation” can be translated into a useful principle: one cannot substantiate causal claims from associations alone, even at the population level—behind every causal conclusion there must lie some causal assumption that is not testable in observational studies. (Pearl 2003) (por si...) - estandarizar el riesgo de la distribución conjunta (joint distribution) de los confusores en todo el seguimiento densidad conjunta de los datos observados para generar resultados potenciales bajo distintos escenarios de resultados --- ## Causalidad (2) **SUTVA** `\(^{[9; 10; 11]}\)` <img src="data:image/png;base64,#./_figs/sutva.jpg" width="300" style="display: block; margin: auto auto auto 0;" /> -- - *Consistencia* -- - *Intercambiabilidad* -- - *Positividad* ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** (2) *#_#_#_#_#_#_#_#_#_#_ potential outcomes framework (Neyman, 1923; Rubin, 1974). ATE = E[Y (1) -Y(0)] - RCTs se puede con la aleatorización. EN estudios observacionales es aún más importante inspeccionar en qué medida se cumplen estos supuestos. stable-unit-treatment-value assumption (SUTVA)” (Rubin 1980). La validez de nuestras inferencias causales requieren dichas condiciones: - Nos obliga a pensar el mecanismo de asignación (assignment mechanism) para definir efectos causales para identificarlos o Consistencia: Si un individuo tiene una exposición A=a, entonces su resultado Y = Y(a). ¿Todos recibieron el tratamiento en medida similar? ¿Lo mismo ocurre con quienes no lo recibieron?  Nos permite homologar los resultados de quienes recibieron una exposición al resultado potencial que se habría observado con ese valor de exposición o Positividad: ¿La distribución de las covariables L en tratados es similar con no-tratados?, la probabilidad de asignación está entre 0 y 1, por lo que se puede tener tratados y no tratados. Cuando hay expuestos y no expuestos en todos los niveles de confusión  Esto a menudo nos e ve en las regresiones. o Intercambiabilidad: Debiese ser la primera, pero la dejé al final por ser la más importante. Fuera de contextos aleatorios, la intercambiabilidad será condicional a los confusores ¿Son comparables los grupos, dada la confusión?  ¿Los tratados y no tratados habrían experimentado en promedio el mismo riesgo de muerte de haber recibido el mismo nivel de tratamiento?  ¿Sigue siendo informativo para el resultado observado saber si son tratados o no, condicional a los niveles de confusión?  los resultados son independientes de la asignación a tratamiento dado que se tienen las covariables - Sesgos como el error de medición, variables no-observadas también podrían afectar el cumplimiento de estos supuestos. --- ## Ajuste confusión - **Estratificación / Regresiones** - **Estandarización** `\(^{[12; 5; 13]}\)` ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ **NO CREO NECESARIO GRAFICAR LA DIFERENCIA ENTRE UNO Y OTRO (AL FINAL FUE MUY DIFÍCIL, PORQUE LOS EJEMPLOS NO ERAN SUFICIENTEMENTE BUENOS** - Asumiendo en gran parte SUTVA, abordaremos las estrategias analíticas para cerrar puertas traseras que pudiesen confundir el efecto del tratamiento en resultado. - La manera que los manejamos los (Lit) puede introducir sesgo (2021 04 21 Marginal Structural Models) - Están los modelos condicionales (Vansteelandt2011) y los basados en la estandarización. - Los primeros agrupan o Estratificación: El riesgo causal específico al estrato es igual al asociacional. Por ejemplo, las regresiones hacen un análisis por cada estrato y de ahí se agrupan; las medidas de resumen (centrar en los niveles de la covariable), o estratificar (restringir el análisis a un nivel de una covariable) (Hernan CURSO) o Estandarización: pseudio-poblacion en que la exposición es independiente de confusores, permitiendo la estimación de parámetros marginales estructurales, es decir, en toda la población de haber recibido el tratamiento. - No obstante, cuando se quiere analizar efectos longitudinales. No basta con el cumplimiento de los supuestos planteados arriba. Muchas exposiciones en epidemiología y salud pública son tiempo-variables (Mansournia2017) - Se suman otras fuentes de confusión y sesgo. Es aquí que son de utilidad los modelos marginales estructurales. Por ejemplo, se generan sesgos de selección en el seguimiento, se introducen confusores tiempo-dependientes (valores que toman más de un valor en el tiempo, algunos producto del tratamiento anterior - - treatment confounder feedback) o tratamientos que son dependientes del tiempo, ajustes que van en un camino desde la exposición al resultado (condicioanr en mediadores) (sobreajuste) u estratificación por colliders - De ahí que la intercambiabilidad pasa a ser secuencial (no sólo en el vector de confusión) Ya0,a1 _||_ A0 --- ## Marginal structural models - ¿Por qué el nombre? - Tres principales: - IPTW - Fórmula G paramétrica - Estimación G ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - IPTW (probabilidad inversa de asignación a tratamiento). Es un modelo en 2 etapas, en que primero estimamos la probabilidad de recibir tratamiento en base a un vector de confusores L. Luego en una segunda etapa se estima el modelo de regresión como tal, pero ya no existiendo la asociación entre confusor y tratamiento, ya que la pseudopoblación ponderada se encuentra balanceada en estos términos. Lo problemático de este modelo es que algunas veces los tratados y no tratados tienen características muy distintas en los niveles de L, lo que puede generar problemas para asignar pesos (ej., asignaría pesos muy extremos e irrealistas). - Fórmula G paramétrica: Tiene un sesgo llamado G-null paradox, en que tiende a rechazar H0 en muestras muy grandes. - La fórmula G es intensiva computacionalmente, al igual la estimación G, requieren más programación y e --- ## Aplicación (1) Se generó una base de datos de 1.000 observaciones con las siguientes características `\(^{[14]}\)` : - t1 = `Tratamiento residencial Tiempo 0 ` ( `\(A_0\)` ) ` ∼ Bernoulli(p=.23) - t2 = `Tratamiento residencial Tiempo 1 ` ( `\(A_1\)` ) ` ∼ Bernoulli(p=.22) - t3 = `Tratamiento residencial Tiempo 2 ` ( `\(A_2\)` ) ` ∼ Bernoulli(p=.15) - c1 = `Policonsumo` ∼ Bernoulli(p=.32) - c2 = `Edad de inicio consumo de sustancias ` ∼ N(3,0.5) - dt1 = `Meses sin readmisión Primer tto.` ∼ Poisson(λ=4) - dt2 = `Meses sin readmisión Segundo tto.` ∼ Poisson(λ=3) - dt3 = `Meses sin readmisión Tercer tto.` ∼ Poisson(λ=2) - v1 = `Confusor tiempo-dependiente (t1)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,basal} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.99})\)` - v2 = `Confusor tiempo-dependiente (t2)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,(2^{do})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.55})\)` - v3 = `Confusor tiempo-dependiente (t3)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,(3^{er})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.33})\)` ` .details-code[ ```r #Semilla para replicar los resultados set.seed(2125) #definir confusor c1 (policonsumo al ingreso==1) policonsumo <-rbinom(1000,1,0.32) #definir confusor c2 (edad de inicio de consumo) edad_ini <-round(exp(rnorm(1000, 3, 0.5))) #definir tratamiento residencial al tiempo 0 t1 <-rbinom(1000,1,0.23) #definir tratamiento residencial al tiempo 1 t2 <-rbinom(1000,1,0.22) #definir tratamiento residencial al tiempo 2 t3 <-rbinom(1000,1,0.15) #meses libre de readmisión T1 dt1 <-rpois(1000,4) #meses libre de readmisión T2 dt2 <-rpois(1000,3) #meses libre de readmisión T3 dt3 <-rpois(1000,1) #definir el confusor tiempo-dependiente v1 como una función de t1 y dt1 v1 <- (0.4*t1 + 0.78*dt1 + rnorm(1000, 0, sqrt(0.99))) + 5 #definir el confusor tiempo-dependiente v2 como una función de t2 y dt2 v2 <- (0.4*t2 + 0.78*dt2 + rnorm(1e3, 0, sqrt(0.55))) + 5 #definir el confusor tiempo-dependiente v3 como una función de t3 y dt3 v3 <- (0.4*t3 + 0.78*dt3 + rnorm(1e3, 0, sqrt(0.33))) + 5 datos <- data.frame(policonsumo, edad_ini, v1, v2, v3, t1, t2, t3, dt1, dt2, dt3) %>% # Generamos una columna id para cada sujeto dplyr::mutate(id = row_number()) %>% # Convertimos la base de datos desde formato ancho a largo de todas las variables tiempo-dependientes. Aclaramos que la estructura de los nombres lleva inscrito primero el valor de interés, seguido por el tiempo/ola de aplicación. Y obtenemos el tiempo tomando el valor dígito de cada columna tidyr::pivot_longer(cols=c("v1", "v2", "v3", "t1", "t2", "t3", "dt1", "dt2", "dt3"), names_to = c('.value',"time"), names_pattern = "(.*)(\\d+)$") %>% # Ordenamos la base de datos por id y tiempo/ola dplyr::arrange(id, time) %>% # reordenamos las columnas dplyr::select(id, time, policonsumo, edad_ini, t, v, dt) %>% # formateamos la base de datos como data.frame as.data.frame() %>% dplyr::mutate(time=as.numeric(time)) #Primeros 4 filas head(datos,4) %>% knitr::kable("html", caption="Muestra de la base de datos (Primeras 4 filas)") %>% kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 10) ``` <table class=" lightable-classic" style='font-size: 10px; font-family: "Arial Narrow", "Source Sans Pro", sans-serif; margin-left: auto; margin-right: auto;'> <caption style="font-size: initial !important;">Muestra de la base de datos (Primeras 4 filas)</caption> <thead> <tr> <th style="text-align:right;"> id </th> <th style="text-align:right;"> time </th> <th style="text-align:right;"> policonsumo </th> <th style="text-align:right;"> edad_ini </th> <th style="text-align:right;"> t </th> <th style="text-align:right;"> v </th> <th style="text-align:right;"> dt </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 10.633698 </td> <td style="text-align:right;"> 7 </td> </tr> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 6.549493 </td> <td style="text-align:right;"> 3 </td> </tr> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 7.595415 </td> <td style="text-align:right;"> 3 </td> </tr> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 44 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 9.276188 </td> <td style="text-align:right;"> 2 </td> </tr> </tbody> </table> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Se asume que todos cursaron el tratamiento y no hay censura, que la velocidad es constante y los riesgos son constantes y proporcionales (Pierce) Auto-regressive (AR1) correlation structure ... correlation to decay as the outcome values were farther away from the time of interest. --- ## Aplicación (2) .details-code[ ```r #```{r, echo=F, dev.args = list(bg = 'transparent'), fig.align="center", out.width=550, error=T} #fig.width = 7, fig.height = 5 muy grande dag34 <- dagitty('dag { bb="0,0,1,1" L1_policonsumo [pos="0.117,0.205"] L2_edad_ini [pos="0.222,0.112"] dt1 [outcome,pos="0.294,0.491"] dt2 [outcome,pos="0.572,0.491"] dt3 [outcome,pos="0.838,0.491"] t1 [exposure,pos="0.092,0.491"] t2 [exposure,pos="0.380,0.491"] t3 [exposure,pos="0.657,0.491"] v1 [pos="0.289,0.686"] v2 [pos="0.576,0.668"] v3 [pos="0.849,0.657"] L1_policonsumo -> dt1 L1_policonsumo -> dt2 L1_policonsumo -> dt3 L1_policonsumo -> t1 L1_policonsumo -> t2 L1_policonsumo -> t3 L2_edad_ini -> dt1 L2_edad_ini -> dt2 L2_edad_ini -> dt3 L2_edad_ini -> t1 L2_edad_ini -> t2 L2_edad_ini -> t3 t1 -> dt1 t1 -> t2 [pos="0.210,0.652"] t2 -> dt2 t2 -> t3 [pos="0.514,0.677"] t3 -> dt3 v1 -> dt1 v1 -> t1 v1 -> t2 v2 -> dt2 v2 -> t2 v2 -> t3 v3 -> dt3 v3 -> t3 dt1 -> dt2 dt2 -> dt3 }') tidy_dag34 <- tidy_dagitty(dag34) %>% dplyr::mutate(label=dplyr::case_when(grepl("A0",as.character(name))~"A0", T~as.character(name))) %>% dplyr::mutate(label2=dplyr::case_when(grepl("L|v",name)~"adj",grepl("LM",name)~"white", T~"black")) %>% dplyr::mutate(adjusted=factor(dplyr::case_when(grepl("L|v",name)~"adjusted",T~"unadjusted"))) edge_function <- ggdag:::edge_type_switch("link_arc") dag34_plot<- ggdag:::if_not_tidy_daggity(tidy_dag34) %>% ggdag:::node_status() %>% ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend, yend = yend, color = status, shape=factor(adjusted)))+ #edge_function()+ scale_adjusted()+ ggdag:::breaks(c("exposure", "outcome","latent"))+ theme_dag()+ geom_dag_edges_arc(curvature = c(rep(0,12),5,5,0,5,0,5,rep(0,9)))+ #14 y 16 de 24 ggdag:::geom_dag_point(size = 16)+ ggdag:::geom_dag_label_repel(ggplot2::aes_string(label = "label", fill = "status"), size = 4.88, col = "white", show.legend = FALSE)+ scale_shape_manual(values = c(15, 16), name="Ajustado", labels=c("Sí", "No"))+ scale_fill_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+ scale_color_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+#E6E6E6 guides(linetype="none", edge_alpha="none", shape="none")+ guides(color=guide_legend(override.aes = list(arrow = NULL)))+#,guide_colourbar(order = 1) theme(plot.caption = element_text(hjust = 0))+ theme(legend.position = "bottom", aspect.ratio=6/10)+ labs(caption="Nota. Ak= Modalidad (Residencial/Ambulatoria); dtk= Meses libre de readmisión;\nvk= Confusor tiempo-dependiente") dag34_plot ``` <img src="data:image/png;base64,#seminario_files/figure-html/tab-comp0_dag-1.png" height="450" style="display: block; margin: auto;" /> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - E[Treatment=1, time=1 | Treatment=1, time=(time-1)] - Generamos una muestra de 1000 casos con las siguientes características. En este ejemplo hipotético tenemos el proceso de generación de los datos y **queremos estimar el efecto de tratamiento residencial en los meses libre de readmisión** - Lo anterior, aplicando un modelo de ponderación inversa de la probabilidad de asignación de tratamiento Estimamos el modelo de probabilidad inversa con el numerador (elementos invariantes en el tiempo) y el denominador todos los confusores de interés. , y luego un modelo marginal estructural en base a un modelo - GEE (generalized estimating ecuations) con errores estándar robustos poisson, con estructura de correlación autorregresiva 1 para datos anidados por paciente (para reflejar la dependencia en los resultados del tratamiento anterior). --- ## Aplicación (3) . - Se estima una ponderación inversa al tratamiento que tiene en cuenta exposición y/o confusores tiempo-dependientes `\(^{[15]}\)` .details-code[ ```r #Paquetes necesarios para generar los análisis library(geepack) library(survey) library(ipw) library(reshape) library(MuMIn) #Estimamos el modelo de probabilidad inversa con el numerador (elementos invariantes en el tiempo) y el denominador todos los confusores de interés. w <- ipwtm( exposure = t, family = "binomial", link = "logit", # Time invariant stuff numerator = ~ factor(policonsumo) + edad_ini, # All confounders denominator = ~ v + factor(policonsumo) + edad_ini, id = id, timevar=time, type="first", data = datos) summary(w$ipw.weights) #Se incorporan los pesos a la base de datos datos<- dplyr::bind_cols(datos, ipw=w$ipw.weights) ########################################################################################### # Para obtener las salidas de los coeficientes ########################################################################################### summary.geem <- function(object, ...) { Coefs <- matrix(NA,nrow=length(object$beta),ncol=5) Coefs[,1] <- c(object$beta) naive <- is.character(object$var) if(!naive && any(diag(object$var) < 0) ){ naive <- TRUE warning("Some elements of robust variance estimate < 0. Reporting model based SE.") } Coefs[,2] <- sqrt(diag(object$naiv.var)) if(naive){Coefs[,3] <- rep(NA, length(object$beta))}else{Coefs[,3] <- sqrt(diag(object$var))} if(naive){Coefs[,4] <- Coefs[,1]/Coefs[,2]}else{Coefs[,4] <- Coefs[,1]/Coefs[,3]} Coefs[,5] <- round(2*pnorm(abs(Coefs[,4]), lower.tail=F), digits=8) colnames(Coefs) <- c("Estimates","Model SE","Robust SE", "wald", "p") summ <<- list(beta = Coefs[,1], se.model = Coefs[,2], se.robust = Coefs[,3], wald.test = Coefs[,4], p = Coefs[,5], alpha = object$alpha, corr = object$corr, phi = object$phi, niter = object$niter, clusz = object$clusz, coefnames = object$coefnames, weights=object$weights, biggest.R.alpha = object$biggest.R.alpha) class(summ) <- 'summary.geem' return(summ) } ########################################################################################### # Primer modelo GEE - GEE con errores estándar robustos (no IPTW) ########################################################################################### mod1<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id ,data = data.frame(datos), family=poisson, corstr="ar1") summary(mod1) summ_mod1<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) invisible("Quasipoisson") mod2<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, family = MASS::negative.binomial(2), corstr = "ar1") #The known value of the additional parameter, theta. mod3<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, family = MASS::negative.binomial(1), corstr = "ar1") mod4<-geepack::geese(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, corstr = "ar1", family=poisson) ``` ] - Una vez que se genera el modelo de ecuaciones de estimación generalizadas (GEE) sin considerar la probabilidad inversa de asignación a tratamiento, se contrasta con el modelo en el que se pondera por dicha probabilidad inversa. Se seleccionó el modelo con menor QIC. .details-code[ ```r #https://gist.github.com/avallecam/56af06f46e5544c3af0f46344df20989 #https://stackoverflow.com/questions/49089476/any-updates-to-model-negative-binomial-distribution-data-with-gee-in-r #https://www.rdocumentation.org/packages/geepack/versions/1.3.4/topics/geese #https://stackoverflow.com/questions/13946540/negative-binomial-in-gee?rq=1 mod1b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id ,weights=ipw, data = data.frame(datos), family=poisson, corstr="ar1") summary.geem(mod1b) summ_mod1b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) mod2b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, weights=ipw, data = datos, family = MASS::negative.binomial(2), corstr = "ar1") #The known value of the additional parameter, theta. summary(mod2b) summ_mod2b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) mod3b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, weights=ipw, data = datos, family = MASS::negative.binomial(1), corstr = "ar1") summary(mod3b) summ_mod3b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) #You can use the quasi-likelihood under the independence model criterion (QIC; Pan, 2001) to assess the model fit. model.sel(mod1b, mod2b, mod3b, rank = QIC) ``` ] Por último, una tabla comparativa de las estimaciones del tratamiento. .details-code[ ```r mt<- rbind.data.frame( summ_mod1[which(summ_mod1$cc=="t"),], summ_mod1b[which(summ_mod1$cc=="t"),] ) rownames(mt) <- c("Modelo sin IPTW", "Modelo con IPTW") mt %>% knitr::kable("html", caption="Comparación resultados") %>% kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 10) ``` <table class=" lightable-classic" style='font-size: 10px; font-family: "Arial Narrow", "Source Sans Pro", sans-serif; margin-left: auto; margin-right: auto;'> <caption style="font-size: initial !important;">Comparación resultados</caption> <thead> <tr> <th style="text-align:left;"> </th> <th style="text-align:left;"> cc </th> <th style="text-align:right;"> Estimate </th> <th style="text-align:right;"> Std.err </th> <th style="text-align:right;"> lwr </th> <th style="text-align:right;"> upr </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> Modelo sin IPTW </td> <td style="text-align:left;"> t </td> <td style="text-align:right;"> 0.99 </td> <td style="text-align:right;"> 1.03 </td> <td style="text-align:right;"> 0.94 </td> <td style="text-align:right;"> 1.05 </td> </tr> <tr> <td style="text-align:left;"> Modelo con IPTW </td> <td style="text-align:left;"> t </td> <td style="text-align:right;"> 0.88 </td> <td style="text-align:right;"> 1.03 </td> <td style="text-align:right;"> 0.83 </td> <td style="text-align:right;"> 0.94 </td> </tr> </tbody> </table> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - applying inverse probability of treatment weights (IPTW) to the observations, which are then fitted to a marginal structural model (MSM). - Se generaron los ponderadores de asignación a tratamiento, mediante un modelo de probabilidad inversa de tratamiento con una estructura, de exposición y confusión tiempo-dependientes. Within each unit under observation this function computes inverse probability weights at each time point during follow-up. T - Se probaron las distribuciones binomiales negativa para ver el grado de sobredispersión de los datos, pero “oh!, qué sorpresa”, parece que no existía tal (bueno, yo construí los datos). - El modelo generalizado simple señala que no hay diferencias en el número de meses libre de de reingresar a tratamiento, entre los tratamientos residenciales y ambulatorios ya que los intervalos de confianza se superponen el valor nulo. No obstante, el modelo que presenta ponderación por el tratamiento indica que los tratamientos residenciales presentan 1/0.88 o un **riesgo 14% menor de estar más meses libre de readmisión, por lo que serían readmitidos 14%más rápido.** *#_#_#_#_#_#_#_#_#_#_ (por si...) We simulate a simple dataset with a time-varying exposure. The data, collected at two time points, t = (0, 1), include a binary outcome measured at the end of follow up Y (1), a binary exposure at two time points {X(0),X(1)}, and a binary confounder of the exposure-outcome relationship that is also affected by prior exposure L(1). L(1) is related to the outcome by the common cause U, which we consider to be unmeasured. This particular simulation setup has been used previously to illustrate other causal inference methods.(Robins y Hernan, 2009) It addresses a central feature of the Bayesian g-formula: control of confounding by covariates that may also be affected by exposure. Estimate ipw weights (time-varying) Estimate inverse probability weights to fit marginal structural models, with a time-varying exposure and time-varying confounders. Within each unit under observation this function computes inverse probability weights at each time point during follow-up. The exposure can be binomial, multinomial, ordinal or continuous. Both stabilized and unstabilized weights can be estimated. estimate inverse probability weights (time-varying) using a logistic regression https://search.r-project.org/CRAN/refmans/ipw/html/ipwtm.html --- ## Desafíos - Structural Nested Models - Confusores tiempo-modificados (*time-modified confounding*) `\(^{[16]}\)` - Generar modelos estructurales de análisis de supervivencia (incluir probabilidad inversa de censura) - Variaciones - Modelos doble o triplemente robustos - Machine Learning (bagging, boosting, random forests, neural networks) ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Métodos mixtos: g-formula y IPTW --- class: center, middle # Gracias! <br> <div class="centered"> Contacto: gonzalez.santacruz.andres@gmail.com </div> <br> <br> <br> .down_right[ <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-6-1.png" width="200" style="display: block; margin: auto;" /> ] .down_left[ <img src="data:image/png;base64,#./_style/Logo_nDP_color_hz_en.png" width="200" style="display: block; margin: auto;" /> ] --- ## Referencias . [1] M. A. Hernán, S. Hernández-Díaz, and J. M. Robins. "A structural approach to selection bias.". Eng. In: _Epidemiology (Cambridge, Mass.)_ 15.5 (sept.. 2004), pp. 615-625. ISSN: 1044-3983 (Print). DOI: [10.1097/01.ede.0000135174.63482.43](https://doi.org/10.1097%2F01.ede.0000135174.63482.43). [2] J. Robins. "A new approach to causal inference in mortality studies with a sustained exposure period-application to control of the healthy worker survivor effect". En. In: _"Math. Modelling"_ 7 (1986), pp. 1393-1512. [3] J. Kaufman. _Conceptos de Causalidad y Sesgos en Estudios Observacionales_. 2020. URL: [https://fivethirtyeight.com/features/science-isnt-broken/#part1](https://fivethirtyeight.com/features/science-isnt-broken/#part1) (visited on 07/09/2022). [4] A. Gelman and G. Imbens. _Why ask Why? Forward Causal Inference and Reverse Causal Questions_. Working Paper 19614. National Bureau of Economic Research, nov.. 2013. DOI: [10.3386/w19614](https://doi.org/10.3386%2Fw19614). URL: [http://www.nber.org/papers/w19614](http://www.nber.org/papers/w19614). [5] M. Ã. Hernán and J. M. Robins. _Causal Inference: What If_. Boca Raton: Chapman & Hall/CRC, 2020. URL: [https://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/2019/10/ci_hernanrobins_1oct19.pdf](https://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/2019/10/ci_hernanrobins_1oct19.pdf). [6] J. Pearl. "Statistics and Causal Inference: A Review". In: _TEST: An Official Journal of the Spanish Society of Statistics and Operations Research_ 12 (feb.. 2003), pp. 281-345. DOI: [10.1007/BF02595718](https://doi.org/10.1007%2FBF02595718). --- ## Referencias (2) . [7] A. Muriel, D. Hernández, and V. Abraira. "Modelos estructurales marginales: una herramienta útil que proporciona evidencia a los estudios observacionales". In: _Nefrologia_ 2.7 (2011), pp. 7-13. ISSN: 20137575. DOI: [10.3265/NefrologiaSuplementoExtraordinario.pre2011.Dec.11267](https://doi.org/10.3265%2FNefrologiaSuplementoExtraordinario.pre2011.Dec.11267). URL: [https://www.revistanefrologia.com/es-modelos-estructurales-marginales-una-herramienta-articulo-X2013757511000636](https://www.revistanefrologia.com/es-modelos-estructurales-marginales-una-herramienta-articulo-X2013757511000636). [8] C. Aschwanden. _Science Isnt Broken Its just a hell of a lot harder than we give it credit for_. 2015. URL: [https://fivethirtyeight.com/features/science-isnt-broken/#part1](https://fivethirtyeight.com/features/science-isnt-broken/#part1) (visited on 06/09/2022). [9] D. B. Rubin. "Randomization analysis of experimental data: The Fisher randomization test comment". In: _Journal of the American statistical association_ 75.371 (1980), pp. 591-593. [10] G. W. Imbens and D. B. Rubin. "Causality: The Basic Framework". In: _Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction_. Cambridge University Press, 2015, p. 3–22. DOI: [10.1017/CBO9781139025751.002](https://doi.org/10.1017%2FCBO9781139025751.002). [11] D. Tompsett, S. Vansteelandt, O. Dukes, et al. "gesttools: General Purpose G-Estimation in R". En. In: _Observational Studies_ 8.1 (2022), pp. 1-28. ISSN: 2767-3324. DOI: [10.1353/obs.2022.0003](https://doi.org/10.1353%2Fobs.2022.0003). URL: [https://muse.jhu.edu/article/856403](https://muse.jhu.edu/article/856403) (visited on 07/10/2022). [12] S. Vansteelandt and N. Keiding. "Invited commentary: G-computation-lost in translation?" En. In: _Am J Epidemiol_ 173.7 (mar.. 2011), pp. 739-742. --- ## Referencias (3) . [13] A. I. Naimi, S. R. Cole, and E. H. Kennedy. "An introduction to g methods". In: _International Journal of Epidemiology_ 46.2 (dic.. 2016), pp. 756-762. ISSN: 0300-5771. DOI: [10.1093/ije/dyw323](https://doi.org/10.1093%2Fije%2Fdyw323). eprint: https://academic.oup.com/ije/article-pdf/46/2/756/25421672/dyw323.pdf. URL: [https://doi.org/10.1093/ije/dyw323](https://doi.org/10.1093/ije/dyw323). [14] M. Bounthavong. _Using inverse probability of treatment weights & Marginal structural models to handle time-varying covariates_. 2018. URL: [https://rpubs.com/mbounthavong/IPTW_MSM_Tutorial](https://rpubs.com/mbounthavong/IPTW_MSM_Tutorial) (visited on 06/09/2022). [15] W. M. van der Wal and R. B. Geskus. "ipw: An R Package for Inverse Probability Weighting". In: _Journal of Statistical Software_ 43.13 (2011), pp. 1-23. URL: [http://www.jstatsoft.org/v43/i13/](http://www.jstatsoft.org/v43/i13/). [16] R. W. Platt, E. F. Schisterman, and S. R. Cole. "Time-modified Confounding". In: _American Journal of Epidemiology_ 170.6 (ago.. 2009), pp. 687-694. ISSN: 0002-9262. DOI: [10.1093/aje/kwp175](https://doi.org/10.1093%2Faje%2Fkwp175). eprint: https://academic.oup.com/aje/article-pdf/170/6/687/285709/kwp175.pdf. URL: [https://doi.org/10.1093/aje/kwp175](https://doi.org/10.1093/aje/kwp175).